{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "084f5229",
   "metadata": {},
   "source": [
    "# The Palatable Diet Problem  \n",
    "\n",
    "**Problem Description.** The Diet Problem is the first large-scale optimization problem to be solved with the Simplex algorithm by Jack Laderman in [1947](https://www.mpi-inf.mpg.de/fileadmin/inf/d1/teaching/winter18/Ideen/Materialien/Dantzig-Diet.pdf). The basic formulation of this problem consists of minimizing the cost of a food basket while meeting the specified nutrient requirements. In this notebook, we solve The Palatable Diet Problem (TPDP), where the basic model is extended with a constraint on the food basket palatability. An explicit formula of the palatability constraint is unknown, but we have data on several food baskets and the respective palatability score. First, we define a conceptual model with the *known constraint*. Then, OptiCL is used to learn and embed the palatability constraint.  \n",
    "(*TPDP is part of a larger optimization problem which simultaneously  optimizes  the  food  basket  to  be  delivered,  the  sourcing  plan,  the  delivery  plan, and  the  transfer  modality  of  a  month-long  food  supply in a Wolrd Food Program setting ([Maragno et al., 2021])*).\n",
    "\n",
    "**Objective function:** minimize the total cost of the food basket.  \n",
    "$\\min_{\\boldsymbol{x}} c^\\top \\boldsymbol{x}$\n",
    "\n",
    "*subject to* \n",
    "\n",
    "**Nutritional constraints:** for each nutrient $j\\in\\mathcal{N}$, at least meet the minimum required level.  \n",
    "$ \\sum_{k \\in \\mathcal{K}} nutval_{kj} x_{k} \\geq nutreq_{j}, \\ \\ \\ \\forall l\\in\\mathcal{N},$   \n",
    "**Constraints on sugar and salt.**</font>  \n",
    "$ x_{salt} = 5,$   \n",
    "$ x_{sugar} = 20,$  \n",
    "**Palatability constraints:** the food basket palatability has to be at least equal to $t$.  \n",
    "$ y \\geq t,$  \n",
    "**Learned predictive model:** the palatability is defined using a predictive model.  \n",
    "$ y = \\hat{h}(\\boldsymbol{x}),$   \n",
    "**Non negativity constraints.**  \n",
    "$ x_{k} \\geq 0, \\ \\ \\ \\forall k \\in \\mathcal{K}.$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "eb2ea7b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from imp import reload\n",
    "import numpy as np\n",
    "import math\n",
    "from sklearn.utils.extmath import cartesian\n",
    "import time\n",
    "import sys\n",
    "import os\n",
    "sys.path.append(os.path.abspath('../../src'))\n",
    "import opticl\n",
    "from pyomo import environ\n",
    "from pyomo.environ import *\n",
    "np.random.seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fd816b1",
   "metadata": {},
   "source": [
    "### Data Loading  \n",
    "**nutr_val**: nutritional values for each of the 25 foods  \n",
    "**nutr_req**: 11 nutrition requirements  \n",
    "**cost_p**: vector of procurement costs  \n",
    "**dataset**: dataframe of food basket instances and relative palatability score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "68b75680",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Beans</th>\n",
       "      <th>Bulgur</th>\n",
       "      <th>Cheese</th>\n",
       "      <th>Fish</th>\n",
       "      <th>Meat</th>\n",
       "      <th>CSB</th>\n",
       "      <th>Dates</th>\n",
       "      <th>DSM</th>\n",
       "      <th>Milk</th>\n",
       "      <th>Salt</th>\n",
       "      <th>...</th>\n",
       "      <th>Soya-fortified bulgur wheat</th>\n",
       "      <th>Soya-fortified maize meal</th>\n",
       "      <th>Soya-fortified sorghum grits</th>\n",
       "      <th>Soya-fortified wheat flour</th>\n",
       "      <th>Sugar</th>\n",
       "      <th>Oil</th>\n",
       "      <th>Wheat</th>\n",
       "      <th>Wheat flour</th>\n",
       "      <th>WSB</th>\n",
       "      <th>label</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>398</th>\n",
       "      <td>0.687675</td>\n",
       "      <td>1.257354</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.302104</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.357429</td>\n",
       "      <td>2.823603</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.637964</td>\n",
       "      <td>0.715428</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3833</th>\n",
       "      <td>0.551125</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.117990</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.392274</td>\n",
       "      <td>2.540599</td>\n",
       "      <td>3.414615</td>\n",
       "      <td>0.733300</td>\n",
       "      <td>0.292719</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4836</th>\n",
       "      <td>0.701614</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.094990</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.330808</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.221908</td>\n",
       "      <td>0.336647</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.545864</td>\n",
       "      <td>0.816616</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4572</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>3.832166</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.626751</td>\n",
       "      <td>0.278648</td>\n",
       "      <td>0.132718</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.311117</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.694007</td>\n",
       "      <td>0.794680</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>636</th>\n",
       "      <td>0.039754</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.344293</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.106482</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.160220</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.788790</td>\n",
       "      <td>0.261417</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4931</th>\n",
       "      <td>0.743757</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.115426</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.700000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.479955</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.041521</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.181913</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.210851</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3264</th>\n",
       "      <td>0.779749</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.075277</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.761232</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.392905</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.178567</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.029013</td>\n",
       "      <td>0.293570</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1653</th>\n",
       "      <td>0.668713</td>\n",
       "      <td>3.322546</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.429041</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.450216</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.271382</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.189249</td>\n",
       "      <td>0.813417</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2607</th>\n",
       "      <td>0.460621</td>\n",
       "      <td>1.464405</td>\n",
       "      <td>0.160834</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.757835</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.391194</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.166244</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.382687</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2732</th>\n",
       "      <td>0.929024</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.339896</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.05</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.259267</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.217787</td>\n",
       "      <td>0.647243</td>\n",
       "      <td>0.647392</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5000 rows × 26 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         Beans    Bulgur    Cheese  Fish  Meat       CSB     Dates       DSM  \\\n",
       "398   0.687675  1.257354  0.000000   0.0   0.0  0.000000  0.000000  0.302104   \n",
       "3833  0.551125  0.000000  0.000000   0.0   0.0  0.000000  0.000000  0.117990   \n",
       "4836  0.701614  0.000000  0.000000   0.0   0.0  0.094990  0.000000  0.330808   \n",
       "4572  0.000000  3.832166  0.000000   0.0   0.0  0.000000  0.626751  0.278648   \n",
       "636   0.039754  0.000000  0.344293   0.0   0.0  0.000000  0.000000  0.106482   \n",
       "...        ...       ...       ...   ...   ...       ...       ...       ...   \n",
       "4931  0.743757  0.000000  0.115426   0.0   0.0  0.700000  0.000000  0.479955   \n",
       "3264  0.779749  0.000000  0.075277   0.0   0.0  0.761232  0.000000  0.392905   \n",
       "1653  0.668713  3.322546  0.000000   0.0   0.0  0.429041  0.000000  0.450216   \n",
       "2607  0.460621  1.464405  0.160834   0.0   0.0  0.757835  0.000000  0.391194   \n",
       "2732  0.929024  0.000000  0.000000   0.0   0.0  0.000000  0.000000  0.339896   \n",
       "\n",
       "          Milk  Salt  ...  Soya-fortified bulgur wheat  \\\n",
       "398   0.000000  0.05  ...                          0.0   \n",
       "3833  0.000000  0.05  ...                          0.0   \n",
       "4836  0.000000  0.05  ...                          0.0   \n",
       "4572  0.132718  0.05  ...                          0.0   \n",
       "636   0.000000  0.05  ...                          0.0   \n",
       "...        ...   ...  ...                          ...   \n",
       "4931  0.000000  0.05  ...                          0.0   \n",
       "3264  0.000000  0.05  ...                          0.0   \n",
       "1653  0.000000  0.05  ...                          0.0   \n",
       "2607  0.000000  0.05  ...                          0.0   \n",
       "2732  0.000000  0.05  ...                          0.0   \n",
       "\n",
       "      Soya-fortified maize meal  Soya-fortified sorghum grits  \\\n",
       "398                    0.000000                           0.0   \n",
       "3833                   0.000000                           0.0   \n",
       "4836                   0.000000                           0.0   \n",
       "4572                   0.000000                           0.0   \n",
       "636                    0.000000                           0.0   \n",
       "...                         ...                           ...   \n",
       "4931                   0.041521                           0.0   \n",
       "3264                   0.000000                           0.0   \n",
       "1653                   0.000000                           0.0   \n",
       "2607                   0.000000                           0.0   \n",
       "2732                   0.000000                           0.0   \n",
       "\n",
       "      Soya-fortified wheat flour  Sugar       Oil     Wheat  Wheat flour  \\\n",
       "398                          0.0    0.2  0.357429  2.823603     0.000000   \n",
       "3833                         0.0    0.2  0.392274  2.540599     3.414615   \n",
       "4836                         0.0    0.2  0.221908  0.336647     0.000000   \n",
       "4572                         0.0    0.2  0.311117  0.000000     0.000000   \n",
       "636                          0.0    0.2  0.160220  0.000000     0.000000   \n",
       "...                          ...    ...       ...       ...          ...   \n",
       "4931                         0.0    0.2  0.181913  0.000000     0.000000   \n",
       "3264                         0.0    0.2  0.178567  0.000000     0.000000   \n",
       "1653                         0.0    0.2  0.271382  0.000000     0.000000   \n",
       "2607                         0.0    0.2  0.166244  0.000000     0.000000   \n",
       "2732                         0.0    0.2  0.259267  0.000000     1.217787   \n",
       "\n",
       "           WSB     label  \n",
       "398   0.637964  0.715428  \n",
       "3833  0.733300  0.292719  \n",
       "4836  0.545864  0.816616  \n",
       "4572  0.694007  0.794680  \n",
       "636   0.788790  0.261417  \n",
       "...        ...       ...  \n",
       "4931  0.000000  0.210851  \n",
       "3264  0.029013  0.293570  \n",
       "1653  0.189249  0.813417  \n",
       "2607  0.000000  0.382687  \n",
       "2732  0.647243  0.647392  \n",
       "\n",
       "[5000 rows x 26 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nutr_val = pd.read_excel('processed-data/Syria_instance.xlsx', sheet_name='nutr_val', index_col='Food', engine='openpyxl')\n",
    "nutr_req = pd.read_excel('processed-data/Syria_instance.xlsx', sheet_name='nutr_req', index_col='Type', engine='openpyxl')\n",
    "cost_p = pd.read_excel('processed-data/Syria_instance.xlsx', sheet_name='FoodCost', index_col='Supplier', engine='openpyxl').iloc[0,:]\n",
    "dataset = pd.read_csv('processed-data/WFP_dataset.csv').sample(frac=1)\n",
    "dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd465891",
   "metadata": {},
   "source": [
    "# OptiCL: Optimization with Constraint Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f98821c2",
   "metadata": {},
   "source": [
    "## Step 1: Conceptual Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3e0a923a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_conceptual_model(cost_p):\n",
    "    N = list(nutr_val.index)  # foods\n",
    "    M = nutr_req.columns  # nutrient requirements\n",
    "\n",
    "    model = ConcreteModel('TPDP')\n",
    "\n",
    "    '''\n",
    "    Decision variables\n",
    "    '''\n",
    "    model.x = Var(N, domain=NonNegativeReals)  # variables controlling the food basket\n",
    "\n",
    "    '''\n",
    "    Objective function.\n",
    "    '''\n",
    "    def obj_function(model):\n",
    "        return sum(cost_p[food].item()*model.x[food] for food in N)\n",
    "\n",
    "    model.OBJ = Objective(rule=obj_function, sense=minimize)\n",
    "\n",
    "    '''\n",
    "    Nutrients requirements constraint.\n",
    "    '''\n",
    "    def constraint_rule1(model, req):\n",
    "        return sum(model.x[food] * nutr_val.loc[food, req] for food in N) >= nutr_req[req].item()\n",
    "    model.Constraint1 = Constraint(M, rule=constraint_rule1)\n",
    "    '''\n",
    "    Sugar constraint\n",
    "    '''\n",
    "    def constraint_rule2(model):\n",
    "        return model.x['Sugar'] == 0.2\n",
    "    model.Constraint2 = Constraint(rule=constraint_rule2)\n",
    "    '''\n",
    "    Salt constraint\n",
    "    '''\n",
    "    def constraint_rule3(model):\n",
    "        return model.x['Salt'] == 0.05\n",
    "    model.Constraint3 = Constraint(rule=constraint_rule3)\n",
    "    \n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32f273e3",
   "metadata": {},
   "source": [
    "## Step 2: Data Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "10151c52",
   "metadata": {},
   "outputs": [],
   "source": [
    "y = dataset['label']\n",
    "X = dataset.drop(['label'], axis=1, inplace=False)\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.8, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d0ba744",
   "metadata": {},
   "source": [
    "## Step 3: Learn the predictive models\n",
    "'alg_list' specifies the list of algorithms that you will consider in the training pipeline. If you have the InterpretableAI license, you can include **iai** (Optimal Trees with Hyperplanes) or **iai-single** (Optimal Trees with single feature splits) in the list. If using IAI, you must specify the metric as 'r2'. Otherwise, the default metric is 'neg_squared_mse'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "48f1f718",
   "metadata": {},
   "outputs": [],
   "source": [
    "version = 'TPDP_v1'\n",
    "alg_list = ['mlp', 'linear', 'gbm', 'cart', 'rf', 'svm']\n",
    "outcome_list = ['palatability']  # Constraint to be learned"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e884eb74",
   "metadata": {},
   "source": [
    "#### Train models (or skip if pre-saved)  \n",
    "The training will use only regression models. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "80fe98fc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Learning a constraint for palatability\n",
      "Training mlp\n",
      "------------- Initialize grid  ----------------\n",
      "------------- Running model  ----------------\n",
      "Algorithm = mlp, metric = None\n",
      "saving... results/mlp_palatability_trained.pkl\n",
      "------------- Model evaluation  ----------------\n",
      "-------------------training evaluation-----------------------\n",
      "Train MSE: 0.011064680927472394\n",
      "Train R2: 0.7832510216353195\n",
      "-------------------testing evaluation-----------------------\n",
      "Test MSE: 0.018394522665774987\n",
      "Test R2: 0.6418989069753125\n",
      "\n",
      "Training linear\n",
      "------------- Initialize grid  ----------------\n",
      "------------- Running model  ----------------\n",
      "Algorithm = linear, metric = None\n",
      "saving... results/linear_palatability_trained.pkl\n",
      "------------- Model evaluation  ----------------\n",
      "-------------------training evaluation-----------------------\n",
      "Train MSE: 0.04494230882369555\n",
      "Train R2: 0.11961315588417554\n",
      "-------------------testing evaluation-----------------------\n",
      "Test MSE: 0.04647042488484837\n",
      "Test R2: 0.09532254536028606\n",
      "\n",
      "Training gbm\n",
      "------------- Initialize grid  ----------------\n",
      "------------- Running model  ----------------\n",
      "Algorithm = gbm, metric = None\n",
      "saving... results/gbm_palatability_trained.pkl\n",
      "------------- Model evaluation  ----------------\n",
      "-------------------training evaluation-----------------------\n",
      "Train MSE: 0.0031460526547302417\n",
      "Train R2: 0.9383711375624779\n",
      "-------------------testing evaluation-----------------------\n",
      "Test MSE: 0.008673420002089097\n",
      "Test R2: 0.8311474975760436\n",
      "\n",
      "Training cart\n",
      "------------- Initialize grid  ----------------\n",
      "------------- Running model  ----------------\n",
      "Algorithm = cart, metric = None\n",
      "saving... results/cart_palatability_trained.pkl\n",
      "------------- Model evaluation  ----------------\n",
      "-------------------training evaluation-----------------------\n",
      "Train MSE: 0.010676479379551228\n",
      "Train R2: 0.7908556050356943\n",
      "-------------------testing evaluation-----------------------\n",
      "Test MSE: 0.015873275949482357\n",
      "Test R2: 0.6909820618521263\n",
      "\n",
      "Training rf\n",
      "------------- Initialize grid  ----------------\n",
      "------------- Running model  ----------------\n",
      "Algorithm = rf_shallow, metric = None\n",
      "saving... results/rf_shallow_palatability_trained.pkl\n",
      "------------- Model evaluation  ----------------\n",
      "-------------------training evaluation-----------------------\n",
      "Train MSE: 0.013386400050741746\n",
      "Train R2: 0.7377702480534164\n",
      "-------------------testing evaluation-----------------------\n",
      "Test MSE: 0.016551194495138304\n",
      "Test R2: 0.6777844716459509\n",
      "\n",
      "Training svm\n",
      "------------- Initialize grid  ----------------\n",
      "------------- Running model  ----------------\n",
      "Algorithm = svm, metric = None\n",
      "saving... results/svm_palatability_trained.pkl\n",
      "------------- Model evaluation  ----------------\n",
      "-------------------training evaluation-----------------------\n",
      "Train MSE: 0.01814469418771081\n",
      "Train R2: 0.6445587582954106\n",
      "-------------------testing evaluation-----------------------\n",
      "Test MSE: 0.019586693402097288\n",
      "Test R2: 0.6186899522497088\n",
      "\n",
      "Saving the performance...\n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "task_type = 'continuous'\n",
    "performance = pd.DataFrame()\n",
    "\n",
    "if not os.path.exists('results/'):\n",
    "    os.makedirs('results/')\n",
    "\n",
    "for outcome in outcome_list:\n",
    "    print(f'Learning a constraint for {outcome}')\n",
    "\n",
    "    for alg in alg_list:\n",
    "        if not os.path.exists('results/%s/' % alg):\n",
    "            os.makedirs('results/%s/' % alg)\n",
    "        print(f'Training {alg}')\n",
    "        s = 1\n",
    "\n",
    "        ## Run shallow/small version of RF\n",
    "        alg_run = 'rf_shallow' if alg == 'rf' else alg\n",
    "\n",
    "        m, perf = opticl.run_model(X_train, y_train, X_test, y_test, alg_run, outcome, task = task_type,\n",
    "                               seed = s, cv_folds = 5, \n",
    "                               # metric = 'r2',\n",
    "                               save = False\n",
    "                              )\n",
    "\n",
    "        ## Save model\n",
    "        constraintL = opticl.ConstraintLearning(X_train, y_train, m, alg)\n",
    "        constraint_add = constraintL.constraint_extrapolation(task_type)\n",
    "        constraint_add.to_csv('results/%s/%s_%s_model.csv' % (alg, version, outcome), index = False)\n",
    "\n",
    "        ## Extract performance metrics\n",
    "        try:\n",
    "            perf['auc_train'] = roc_auc_score(y_train >= threshold, m.predict(X_train))\n",
    "            perf['auc_test'] = roc_auc_score(y_test >= threshold, m.predict(X_test))\n",
    "        except: \n",
    "            perf['auc_train'] = np.nan\n",
    "            perf['auc_test'] = np.nan\n",
    "\n",
    "        perf['seed'] = s\n",
    "        perf['outcome'] = outcome\n",
    "        perf['alg'] = alg\n",
    "        perf['save_path'] = 'results/%s/%s_%s_model.csv' % (alg, version, outcome)\n",
    "        \n",
    "            \n",
    "        perf.to_csv('results/%s/%s_%s_performance.csv' % (alg, version, outcome), index = False)\n",
    "        \n",
    "        performance = performance.append(perf)\n",
    "        print()\n",
    "print('Saving the performance...')\n",
    "performance.to_csv('results/%s_performance.csv' % version, index = False)\n",
    "print('Done!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b937d3a",
   "metadata": {},
   "source": [
    "## Step 4: Predictive model selection and Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "644653d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "What is the palatability threshold that you want to use in the constraint? The default is 0.5.\n",
      "Choose in the range (0, 1): 0.5\n"
     ]
    }
   ],
   "source": [
    "question2 = 0.5\n",
    "question2 = input('What is the palatability threshold that you want to use in the constraint? The default is 0.5.\\nChoose in the range (0, 1): ')\n",
    "if float(question2) > 1:\n",
    "    question2 = 0.5\n",
    "elif float(question2) < 0:\n",
    "    question2 = 0.5\n",
    "threshold = float(question2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "99e47908",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>save_path</th>\n",
       "      <th>seed</th>\n",
       "      <th>cv_folds</th>\n",
       "      <th>task</th>\n",
       "      <th>parameters</th>\n",
       "      <th>best_params</th>\n",
       "      <th>valid_score</th>\n",
       "      <th>train_score</th>\n",
       "      <th>train_r2</th>\n",
       "      <th>test_score</th>\n",
       "      <th>test_r2</th>\n",
       "      <th>outcome</th>\n",
       "      <th>alg</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>results/mlp/TPDP_v1_palatability_model.csv</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>continuous</td>\n",
       "      <td>{'hidden_layer_sizes': [(10,), (20,), (50,), (...</td>\n",
       "      <td>{'hidden_layer_sizes': (100,)}</td>\n",
       "      <td>-0.019855</td>\n",
       "      <td>0.011065</td>\n",
       "      <td>0.783251</td>\n",
       "      <td>0.018395</td>\n",
       "      <td>0.641899</td>\n",
       "      <td>palatability</td>\n",
       "      <td>mlp</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>results/linear/TPDP_v1_palatability_model.csv</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>continuous</td>\n",
       "      <td>{'alpha': [0.1, 1, 10, 100, 1000], 'l1_ratio':...</td>\n",
       "      <td>{'alpha': 0.1, 'l1_ratio': 0.1}</td>\n",
       "      <td>-0.045594</td>\n",
       "      <td>0.044942</td>\n",
       "      <td>0.119613</td>\n",
       "      <td>0.046470</td>\n",
       "      <td>0.095323</td>\n",
       "      <td>palatability</td>\n",
       "      <td>linear</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>results/gbm/TPDP_v1_palatability_model.csv</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>continuous</td>\n",
       "      <td>{'learning_rate': [0.01, 0.025, 0.05, 0.075, 0...</td>\n",
       "      <td>{'learning_rate': 0.2, 'max_depth': 5, 'n_esti...</td>\n",
       "      <td>-0.008329</td>\n",
       "      <td>0.003146</td>\n",
       "      <td>0.938371</td>\n",
       "      <td>0.008673</td>\n",
       "      <td>0.831147</td>\n",
       "      <td>palatability</td>\n",
       "      <td>gbm</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>results/cart/TPDP_v1_palatability_model.csv</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>continuous</td>\n",
       "      <td>{'max_depth': [3, 4, 5, 6, 7, 8, 9, 10], 'min_...</td>\n",
       "      <td>{'max_depth': 10, 'max_features': 1.0, 'min_sa...</td>\n",
       "      <td>-0.016014</td>\n",
       "      <td>0.010676</td>\n",
       "      <td>0.790856</td>\n",
       "      <td>0.015873</td>\n",
       "      <td>0.690982</td>\n",
       "      <td>palatability</td>\n",
       "      <td>cart</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>results/rf/TPDP_v1_palatability_model.csv</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>continuous</td>\n",
       "      <td>{'n_estimators': [10, 25], 'max_features': ['a...</td>\n",
       "      <td>{'max_depth': 4, 'max_features': 'auto', 'n_es...</td>\n",
       "      <td>-0.016288</td>\n",
       "      <td>0.013386</td>\n",
       "      <td>0.737770</td>\n",
       "      <td>0.016551</td>\n",
       "      <td>0.677784</td>\n",
       "      <td>palatability</td>\n",
       "      <td>rf</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>results/svm/TPDP_v1_palatability_model.csv</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>continuous</td>\n",
       "      <td>{'C': [0.1, 1, 10, 100]}</td>\n",
       "      <td>{'C': 100}</td>\n",
       "      <td>-0.019034</td>\n",
       "      <td>0.018145</td>\n",
       "      <td>0.644559</td>\n",
       "      <td>0.019587</td>\n",
       "      <td>0.618690</td>\n",
       "      <td>palatability</td>\n",
       "      <td>svm</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                       save_path  seed  cv_folds        task  \\\n",
       "0     results/mlp/TPDP_v1_palatability_model.csv     1         5  continuous   \n",
       "1  results/linear/TPDP_v1_palatability_model.csv     1         5  continuous   \n",
       "2     results/gbm/TPDP_v1_palatability_model.csv     1         5  continuous   \n",
       "3    results/cart/TPDP_v1_palatability_model.csv     1         5  continuous   \n",
       "4      results/rf/TPDP_v1_palatability_model.csv     1         5  continuous   \n",
       "5     results/svm/TPDP_v1_palatability_model.csv     1         5  continuous   \n",
       "\n",
       "                                          parameters  \\\n",
       "0  {'hidden_layer_sizes': [(10,), (20,), (50,), (...   \n",
       "1  {'alpha': [0.1, 1, 10, 100, 1000], 'l1_ratio':...   \n",
       "2  {'learning_rate': [0.01, 0.025, 0.05, 0.075, 0...   \n",
       "3  {'max_depth': [3, 4, 5, 6, 7, 8, 9, 10], 'min_...   \n",
       "4  {'n_estimators': [10, 25], 'max_features': ['a...   \n",
       "5                           {'C': [0.1, 1, 10, 100]}   \n",
       "\n",
       "                                         best_params  valid_score  \\\n",
       "0                     {'hidden_layer_sizes': (100,)}    -0.019855   \n",
       "1                    {'alpha': 0.1, 'l1_ratio': 0.1}    -0.045594   \n",
       "2  {'learning_rate': 0.2, 'max_depth': 5, 'n_esti...    -0.008329   \n",
       "3  {'max_depth': 10, 'max_features': 1.0, 'min_sa...    -0.016014   \n",
       "4  {'max_depth': 4, 'max_features': 'auto', 'n_es...    -0.016288   \n",
       "5                                         {'C': 100}    -0.019034   \n",
       "\n",
       "   train_score  train_r2  test_score   test_r2       outcome     alg  \n",
       "0     0.011065  0.783251    0.018395  0.641899  palatability     mlp  \n",
       "1     0.044942  0.119613    0.046470  0.095323  palatability  linear  \n",
       "2     0.003146  0.938371    0.008673  0.831147  palatability     gbm  \n",
       "3     0.010676  0.790856    0.015873  0.690982  palatability    cart  \n",
       "4     0.013386  0.737770    0.016551  0.677784  palatability      rf  \n",
       "5     0.018145  0.644559    0.019587  0.618690  palatability     svm  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "constraints_embed = ['palatability']\n",
    "objectives_embed = {}\n",
    "performance = pd.read_csv('results/%s_performance.csv' % version)\n",
    "performance.dropna(axis='columns')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b00d0715",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "        outcome model_type                                   save_path  \\\n",
      "0  palatability        gbm  results/gbm/TPDP_v1_palatability_model.csv   \n",
      "\n",
      "         task  objective  \n",
      "0  continuous          0  \n"
     ]
    }
   ],
   "source": [
    "model_master = opticl.model_selection(performance, constraints_embed, objectives_embed)\n",
    "\n",
    "model_master['lb'] = threshold\n",
    "model_master['ub'] = None\n",
    "model_master['SCM_counterfactuals'] = None\n",
    "model_master['features'] = [[col for col in X.columns]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "10d32284",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Do you want to use the trust region? True\\False: True\n"
     ]
    }
   ],
   "source": [
    "trust_region = input('Do you want to use the trust region? True\\False: ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6cadcab9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getSolution(model, X):\n",
    "    solution = {}\n",
    "    palatability = 0\n",
    "    count = 0\n",
    "    for v in model.getVars():\n",
    "        if 'x[' in v.varName:\n",
    "            solution[list(X.columns)[count]]=[v.x]\n",
    "            print(v.varName)\n",
    "            count += 1\n",
    "    return solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6434a547",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating constraints for the trust region using 5000 samples.\n",
      "... Trust region defined.\n",
      "Embedding constraints for palatability\n",
      "---------------------Solving the optimization problem---------------------\n",
      "The optimal solution is: \n",
      " {'Beans': '16.03g', 'DSM': '0.02g', 'Milk': '56.78g', 'Salt': '5.0g', 'Lentils': '36.03g', 'Maize': '130.81g', 'Chickpeas': '0.34g', 'Sugar': '20.0g', 'Oil': '24.05g', 'Wheat': '215.27g', 'WSB': '70.02g'}\n",
      "The predicted palatability of the optimal solution is 0.502000550485145\n"
     ]
    }
   ],
   "source": [
    "result = {}\n",
    "conceptual_model= init_conceptual_model(cost_p)\n",
    "MIP_final_model = opticl.optimization_MIP(conceptual_model, conceptual_model.x, model_master, X, tr=trust_region=='True')\n",
    "opt = SolverFactory('gurobi')\n",
    "print('---------------------Solving the optimization problem---------------------')\n",
    "results = opt.solve(MIP_final_model) \n",
    "solution = {}\n",
    "for food in  list(nutr_val.index):\n",
    "    if value(MIP_final_model.x[food])*100 > 0.0000001:\n",
    "        solution[food] = str(np.round(value(MIP_final_model.x[food])*100, 2))+'g'\n",
    "print('The optimal solution is: \\n', solution)\n",
    "print(f\"The predicted palatability of the optimal solution is {value(MIP_final_model.y['palatability'])}\")"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
